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Abstract 

The evolution of a spherical single-mass star cluster is followed in detail up 
to core collapse by numerically solving the orbit-averaged two-dimensional Fokker- 
Planck equation in energy-angular momentum space. Velocity anisotropy is allowed 
in the two-dimensional Fokker-Planck model. Using improved numerical codes, the 
evolution has been followed until the central density increased by a factor of 10^^ 
with high numerical accuracy. The numerical results clearly show self-similar evolu¬ 
tion of the core during the late stages of the core collapse. In the self-similar region 
between the isothermal core and the outer halo, the density profile is characterized 
by a power law p{r) oc and the ratio of the one-dimensional tangential ve¬ 
locity dispersion to the radial one is — 0.92. As the core collapse proceeds, 

the collapse rate ^ = tr(0)dln p{0)/dt tends to the limiting value of ^ = 2.9 x 10“^, 
which is 19% smaller than the value for isotropic clusters. When Plummer’s model 
is chosen as the initial condition, the collapse time is about 17.6 times the initial 
half-mass relaxation time. As the result of strong relaxation in the core, the halo 
becomes to be dominated by radial orbits. The degree of anisotropy monotonically 
increases as the radius increases. In the outer halo, the profiles of the density are 
approximated by p oc This work confirms that the generation of velocity 

anisotropy is an important process in collisional stellar systems. 

Key words: Clusters: globular — Fokker-Planck equation — Numerical meth¬ 
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1 Introduction 


Since Cohn (1979, 1980) invented a numerical scheme to solve the orbit-averaged 
Fokker-Planck (hereafter FP) equation by direct integration, the scheme has been widely 
used as a main tool to study the dynamical evolution of globular clusters. Cohn (1979) hrst 
performed direct numerical integration of the time-dependent two-dimensional (hereafter 
2D) FP equation in energy-angular momentum {E, J) space with calculation of a self- 
consistent potential. This work was done after Cohn and Kulsrud (1978), who studied a 
steady-state stellar distribution around a central massive black hole by direct numerical 
integration of the 2D FP equation. Although Cohn (1979) showed the potential power 
of the direct integration scheme, he had to stop the calculation at a relatively early 
stage of gravothermal core collapse due to numerical error in energy conservation. The 
central density of the cluster had increased by only three orders of magnitude as compared 
with its initial value, when the calculation was stopped. Later Cohn (1980) treated one¬ 
dimensional (hereafter ID) energy space FP equation which was obtained by integrating 
the {E, J) space FP equation over J-space. [This i7-space FP equation was previously 
employed by Henon (1961).] The averaging over J-space means that we assume isotropy 
of the velocity distribution, i.e. spherical symmetry in velocity space. 

One of the advantages of using the ID FP equation is that it saves a lot of computation 
time and computer memory space. Another advantage is that numerical error in energy 
conservation is much reduced (Cohn 1980). This reduction of the error is largely owing 
to the adoption of Chang and Cooper’s (1970) hnite-differencing scheme for the ID FP 
equation. The characteristics of the Chang-Cooper scheme are that it ensures particle 
conservation and non-negativeness of the distribution function, and reproduces the ana¬ 
lytic quasi-equilibrium solution exactly. Due to these characteristics the Chang-Cooper 
scheme achieves high accuracy with a relatively small number of meshes. Cohn (1980) 
reported that the secular energy drift rate was reduced by better than a factor of 100 
with the adoption of the Chang-Cooper scheme. In fact he could follow the core collapse 
until the central density increased by twenty orders of magnitude. 

Because of the above-mentioned advantages of the ID FP equation, since the work 
of Cohn (1980), the ID FP equation has been used in most studies on globular cluster 
evolution while the 2D FP equation has been seldom used. The adoption of the ID E- 
space (isotropic) FP equation seems to be quite reasonable to study the core evolution 
which has been a main subject during the past two decades, because strong relaxation 
will enforce isotropy of the velocity distribution in the core. In fact, Cohn (1985) argued 
that anisotropy does not signihcantly alter the core collapse, although he found that 
anisotropy extended well inside the half-mass radius at late times of the core collapse. 
Our understanding of the evolution of globular clusters has improved very much during 
the past hfteen years, which has been largely owing to the use of the simple ID FP 
equation (cf. Goodman 1993). 

The development of anisotropy in the halo is a natural consequence of cluster evolution 
driven by two-body relaxation (Wooley, Robertson 1956; Henon 1971; Spitzer, Hart 1971b; 
Bettwieser, Spurzem 1986). The relaxation is strong in the core because of its high density 
and the strong relaxation continues to produce high-energy stars. The orbits of the high- 
energy stars extends into the halo. Because the density is very low in the halo, the stars 
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travel the halo almost without experiencing collisions and return to the core. Thus the 
high-energy stars scattered out of the core have very radial orbits on the average. Halo 
stars which are initially not on very radial orbits do not go into the core and they evolve 
little. Therefore radial orbits predominate the halo and velocity anisotropy increases as 
the halo grows. 

The penetration of anisotropy even into the inner region was pointed out in several 
early works (Larson 1970; Spitzer, Shull 1975; Duncan, Shapiro 1982; Cohn 1985). Spitzer 
and Shapiro (1975) argued that it was a natural consequence of the gravothermal core 
collapse. The core becomes hotter as it collapses. Stars on more nearly radial orbits 
more reflects this increased temperature than stars on more nearly circular orbits. Thus 
anisotropy in the core is expected to penetrate into the more inner region as the core 
collapses. 

As described above, the development of velocity anisotropy is expected through the 
whole cluster. So far anisotropy was considered in various simulations of the evolution 
of star clusters, e.g. in calculations of moment equations of the FP equation (Larson 
1970; Bettwieser 1983), and in Monte-Carlo calculations of the FP equation (Henon 1971; 
Spitzer, Shull 1975; Duncan, Shapiro 1982). However, direct 2D FP calculations were 
carried out only by Cohn(1979, 1985). Cohn’s calculations were for isolated, single-mass 
clusters without binaries. Direct 2D FP calculations of more realistic globular cluster 
models which include the effects of binaries, stellar mass spectrum, the galactic tidal 
held, etc. have not been reported until now. Therefore how velocity anisotropy changes 
the evolution of such realistic cluster models is not known well. 

Recently various anisotropic gaseous models and higher-order huid-dynamical mod¬ 
els of star clusters have been developed (Bettwieser, Spurzem 1986; Louis 1990; Louis, 
Spurzem 1991; Spurzem 1991; Giersz, Spurzem 1994). These models are composed of 
moment equations of the FP equation with some closure relations. These models have 
made it possible to calculate in detail the evolution of anisotropic star clusters. Some cal¬ 
culations of the post-collapse evolution and of the evolution of multi-mass clusters have 
been already done (Giersz, Spurzem 1994; Spurzem 1994; Spurzem, Takahashi 1995). 
However, there are some uncertainties in these models. It is unclear what closures we 
should choose, and whether local description of relaxation effects is really valid for stellar 
systems. Therefore the reliability of anisotropic gaseous and higher-order fluid-dynamical 
models should be checked by comparing these results with the results of more fundamen¬ 
tal FP models and/or iV-body models. Some such comparisons have been already done 
(Giersz, Spurzem 1994; Spurzem, Takahashi 1995); the results of the anisotropic gaseous 
models are generally in good agreement with those of the other models. However, com¬ 
parison with the 2D FP models has not been done, because any reliable 2D FP code was 
not available. 

Now is the good time to reconsider the direct 2D FP calculations. The sophisticated 
anisotropic gaseous models are available as described above. The comparison between the 
gaseous and the FP models is very useful for checking the reliability of these models and 
for interpreting the results of these models. On the practical sides of computation, the 
capability of computers has greatly developed since the hrst direct 2D FP calculations 
were carried out by Cohn (1979). At present it is possible to carry out the 2D calculations 
on standard workstations . 

A main obstacle to the 2D FP calculations is to develop computational schemes of high 
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accuracy (especially in the energy conservation). Concerning this problem, a numerical 
scheme which may be good for 2D FP calculations was proposed by Takahashi (1993a). 
He utilized a scheme based on a generalized variational method instead of the Chang- 
Cooper scheme to solve the FP equation. Although a classical variational principle for 
the FP equation does not exist, a generalized variational principle based on the concept 
of a local potential (Glansdorff, Prigogine 1971) does. Inagaki and Lynden-Bell (1990) 
gave the local potentials for the orbit-averaged FP equations. The generalized variational 
method was applied by Takahashi and Inagaki (1992) and Takahashi (1993b) to obtain 
self-similar solutions of the ID FP equation. Takahashi (1993a) solved the time-dependent 
ID FP equation by his method and showed that the method achieved an equal numerical 
accuracy with the accuracy achieved by Cohn’s method (or the Chang-Cooper scheme). 
In Takahashi’s method the FP equation is discretized by the hnite element method based 
on the variational method. Takahashi (1993a) suggested that his method may be useful 
for the 2D FP equation as well, because no special technique applicable only for ID cases 
was used in the method. On the other hand, any direct 2D generalization of the Chang- 
Cooper scheme (which preserves all merits of the Chang-Cooper scheme) has not been 
found as far as the author knows. 

The aim of the present paper is to develop reliable numerical schemes for the 2D orbit- 
averaged FP equation and to consider in detail the evolution of spherical and anisotropic 
star clusters. In other words, we improve Cohn’s (1979) work so that we can understand 
the evolution of anisotropic clusters as well as the evolution of isotropic clusters. In the 
present paper, we consider only the pre-collapse evolution of isolated, single-mass clusters. 
In particular we pay attention to the self-similar structure of the core collapse and the 
development of the strongly anisotropic halo. The effects of binary heating, stellar mass 
spectrum, the galactic tidal held, etc. will be considered in future works. 

In section |], the orbit-averaged FP equation in {E, J)-space is described briehy. In 
section the numerical methods to solve the FP equation are described. The details of the 
discretization schemes are given in [Appendix 1|. The accuracy of the present calculations 


is discussed in section |. In section |, the results of the calculations are presented. In 
section |6|, the results are summarized. 


2 The Orbit-Averaged Fokker-Planck Equation 


We study the evolution of spherical one-component star clusters. We dehne a distri¬ 
bution function /(r, v, t) so that /(r, v, t) d^r is the number of stars at time t within 
the volume element d^r centered at r and within the velocity space element d^v centered 
at V. In a steady-state spherical system, / is a function of only energy E and angular 
momentum J per unit mass. Then the distribution function f{E, J) changes only through 
collisional effects. The time scale of the change is the two-body relaxation time and it is 
much longer than the dynamical time scale td in stellar systems with large N (cf. Binney, 
Tremaine 1987; Spitzer 1987). Thus we can assume that dynamical equilibrium is always 
established when we consider the evolution of / caused by collisions. In this case, the 
evolution of / can be described by the orbit-averaged FP equation in {E, J)-space (Cohn 
1979). 

In numerical calculations, it is more convenient to use the scaled angular momentum 
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R instead of J as a basic variable. R is defined by 


R 


.P 

JHE)' 


(1) 


where Jc{E) is the angular momentum of a circular orbit of energy E. Thus R takes all 
values between 0 and 1, independent of E. {R = 0 corresponds to radial orbits and R = 1 
to circular orbits.) The number density N{E,R) in {E,R)-space is given by 


N{E,R) = 4PP{E,R)P^{E)f{E,R) 

^ A{E,R)f{E,R), (2) 


where P{E, R) is the orbital period. As Cohn (1979) showed, the 2D FP equation under 
the hxed gravitational potential (j){r) can be written in a flux conserving form 
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The expressions for coefficients Deei De, etc. are given in Appendix C of Cohn (1979). 

An isotropized distribution function has been used to derive these expressions for the 
coefficients. The isotropized distribution function f{E,r) is dehned as 


fiE,r) 


1 m 


( 5 ) 


where Rjnax{E,r) = 2p[(j)[r) — E]/J^{E) is the maximum allowed value of R for all the 
orbits of energy E which pass through the radius r. (We adopt the opposite of the 
usual sign convention for the potential 0 throughout this paper; thus E = 0 —n^/2.) The 
isotropization makes the computation of the coefficients much easier. It has been generally 
conceived that the use of the isotropized distribution function does not bring signihcant 
inaccuracies, because the coefficients only depend on moments of f{E,R). However, as 
we see below, f{E, R) becomes to strongly depend on i? for ~ 0 as the halo develops. 
Thus the isotropization may bring signihcant error in such cases. On the other hand, 
since the relaxation occurs virtually only in the core where the distribution function is 
almost entirely isotropic, the error may be negligible. Although careful studies on the 
effect of the isotropization is desirable, we leave such studies to the future. 


3 Methods of Calculation 


The framework of our method is the same as that of Cohn’s (1979) method. Cohn’s 
method is made up of two steps, which are the FP step and the Poisson step. In the 
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FP step the distribution function is advanced by solving the FP equation with the grav¬ 
itational potential being held hxed. In the Poisson step the potential is advanced by 
solving Poisson’s equation with the distribution function being held hxed as a function of 
adiabatic invariants. 

A essential difference between our method and Cohn’s method exists only in a dis¬ 
cretization scheme of the FP equation. Here, therefore, we mainly describe procedures in 
the FP step. 

We discretize the 2D FP equation (^) by either the hnite difference method or the 
hnite element method. Concerning the hnite diherence method, a Chang-Cooper-like 
scheme is applied only to the F^-direction. This scheme improves the accuracy of the total 
energy conservation as compared with a simple centered diherence scheme. Concerning 
the hnite element method, it was found that the test and weight functions implied by the 
generalized variational principle (Inagaki, Lynden-Bell 1990; Takahashi 1993) are good 
for the numerical accuracy. The details of the discretization schemes are described in 
Appendix 1| . After the discretization procedure, the set of linear algebraic equations is 
obtained, and it is solved by an appropriate solver of matrix equations. 

The FP equation is solved in a rectangular domain enclosed by boundary lines, E = 
0(0), E = = 0, and i? = 1, where 0(0) is the central potential and the value 

of i^min is chosen to be close to zero. We impose boundary conditions that Fg = 0 
on boundaries E = 0(0), F min and Fr = 0 on boundaries F = 0,1. The only nontrivial 
condition is Fr = 0 on the F = F^in boundary. We choose this condition so that the total 
mass is conserved. When we set the value of Fmin close enough to zero, what boundary 
condition we choose may not seriously affect the results we are interested in. 

We use variables {X, Y) instead of (F, R) in practical calculations. The variable X{E) 
is defined by 

F 


X{E) = In 


( 6 ) 


_2(>(0) -E„-E 

where Eq is an adjustable parameter. This variable was introduced by Cohn (1979). The 
value of Eq was chosen to be equal to the energy of a circular orbit at the core radius. 
The variable Y (R) is defined by 


Y{R) 


In (1 -|- RjR q) 
ln(l + l/Fo) ’ 


(7) 


where Rq is an adjustable parameter such as 0 < Fq 1. We set Rq = 0.01 in standard 
runs. This variable is introduced in order to give good representation to radial orbits. 
The necessity of introducing this variable is clear when we see figure 10a. We set uni¬ 
form meshes in X and Y. We also need a radial mesh and set a uniform mesh in logr 
between r ^in and rmax- In the standard runs, we chose r ^in = 10“®ro, rmax = lOOro, 
and F min = 0(rmax), where tq a length scale parameter of Plummer’s model chosen as an 
initial distribution. 

The time step Afp between two succeeding Poisson steps is chosen in such a way that 
the fractional increase in the central density during Afp is always about 2.5%. The time 
step Atpp used in the FP integration is set to be Afp/10 in the standard runs. 
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4 Numerical Accuracy 


Numerical error, particularly in energy conservation, was a serious problem in Cohn’s 
(1979) 2D FP calculations as described in section [^. Due to the problem, direct 2D FP 
calculations have been seldom performed since then. Therefore we examine the numerical 
accuracy of our calculations before we discuss the results. 

The numbers of grid points in X, Y, and r are denoted by Nx, Ny, and N^, respec¬ 
tively. We set Nx = 201, Ny = 25, and W = 101 in the standard runs. Test runs were 
carried out with other sets of grid numbers. The calculations started from Plummer’s 
model were continued until the central density increased by about 14 orders of magnitude 
in the standard runs. 

Since we used the boundary conditions of vanishing particle-fluxes as described in 
section ^ the total cluster mass should be conserved. In the FP steps, if we use the hnite 
difference scheme described in appendix 1.1, the total mass is exactly conserved up to 
round-off error because the flux conserving hnite-difference scheme is used. However, the 
mass error is caused by the Poisson steps. During the calculation the fractional change 
in the mass was within 0.1%. If we use the hnite element scheme described in appendix 
1.2, the total mass is not exactly conserved in the FP steps unfortunately. In this case 
the fractional change was within 0.3%. 

As Cohn (1979) showed in his Appendix A.II, the total energy should be conserved, 
except for the change due to the energy hux through the E = i^min boundary, by the 
two-step procedure consisting of a FP step and a following Poisson step. The rate of the 
change of the total (binding) energy £ (due to the energy hux through the boundary) is 
given by 


dS Y Y 


U < (AB)" >. / 


( 8 ) 


where < [AE)'^ >t is the orbit-averaged dihusion coefficient concerning the squared 
change of the energy. The hrst term of the right-hand side vanishes because we im¬ 
pose the boundary condition ^^(i^minj-R) = 0. The second term represents a dihusive 
energy hux through the E = F^min boundary. This term is generally nonzero (positive) 
unless /(i^min, R) = 0, even if the particle hux Ee is zero at the E = i^min boundary. 
We might think that this energy change could be negligible since the value of / was very 
small at the E = i^min boundary. This is true in ID i7-space FP calculations. In the 
present 2D calculations, however, strong anisotropy develops in the halo and the value of 
/ is rather large near i? = 0 even at the E = F^min boundary (see hgure 10a). The amount 
of the cumulative energy changes given by equation (||) was about 1% of the total energy 
at the end of the present calculations. Considering these changes, we estimated that the 
numerical error in the total energy were within 1% during the calculations. In early stages 
of the calculations, the rate of the energy change (|]) increased rapidly due to the rapid 
development of anisotropy. In these stages, the error in the energy mainly came from the 
FP steps. The error seemed to be smaller for the hnite element code than for the hnite 
diherence code. In middle stages of the calculations, the energy error was smaller for the 
hnite diherence code when Nx was relatively small (e.g. iVx=151). In late stages of the 
calculations (then only the core actually evolved), the energy error mainly came from the 
Poisson steps and the total (binding) energy steadily decreased. We could not dehnitely 
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conclude which of the two codes was better from a viewpoint of the energy conservation. 
When the simple centered hnite-difference scheme was used, the energy error amounted 
to about 8% at the end of the calculation. 

The reliability of the calculations can be indirectly tested by comparing the results 
obtained by various numerical methods. In the present paper, the two different dis¬ 
cretization schemes were used. The results obtained by the two schemes were compared 
in various points as described in section |^, and they were generally in good agreement. 
This fact supports the reliability of our integration schemes of the FP equation. Only 
noticeable differences between the two schemes appeared in the evolution of quantities ^ 
and 7 dehned by equations ( 0 ) and ([T3|) . These quantities are dehned in the core and 
tends to be constant as the core collapse proceeds (see hgures 4 and 5). When the hnite 
element scheme was used, the limiting constant values of ^ and 7 decreased as the number 
of X(energy)-grids Nx increased: x 10^, 7 ) = (4.06, 0.111), (3.53, 0.109), (3.28, 0.107) 

for Ab(:=101, 151, 201 at the end of each run. When the hnite diherence scheme was used, 
the limiting values varied little as Nx increased, i.e. the results seemed to have already 
converged: (^ x 10^, 7 ) = (2.98, 0.106), (2.95, 0.105), (2.94, 0.104) for Arx=101, 151, 201 
at the end of each run. The limiting values by the hnite element code seem to converge 
to the limiting values by the hnite diherence code as Nx increases. We may regard the 
limiting values of ^ and 7 obtained by the hnite diherence code as the real values. 


5 Results 


Calculations were carried out by using both the hnite diherence and hnite element 
codes. The results by the two (partially) diherent codes were in good agreement except 
on the evolution of ^ and 7 which is described in section Figures shown in this section 
were actually drawn from the results of the calculation by the hnite diherence code. 

The initial condition of our calculations is Plummer’s model where the velocity distri¬ 
bution is isotropic everywhere. We use the standard units of G = M = 1 and £i = 1/4, 
where G is the gravitational, M is the initial total mass, and E\ is the initial total (binding) 
energy. No heat sources are included in the present calculations. 

Figure 1 shows the evolution of the density prohle. The evolving density prohle is 
plotted every 200 potential-recalculation time steps. The central density continues to in¬ 
crease and the core radius continues to decrease. This is due to the well-known gravother- 
mal instability (Antonov 1962; Lynden-Bell, Wood 1968). The self-similar nature of the 
gravothermal core collapse (Lynden-Bell, Eggleton 1980) is clearly seen in hgure 1. 

The radial prohle of the logarithmic density gradient —a = din p/d In r is shown in 
hgure 2. From this hgure we hud that the inner power-law region of p oc develops 

as the core collapse proceeds. This value of a = 2.23 coincides with that found in a ID 
(isotropic) FP calculation (Cohn 1980). Heggie and Stevenson (1988) also found the same 
value of a = 2.23 for the pre-collapse self-similar solution of the isotropic FP equation. 
A 2D (anisotropic) FP calculation by Cohn (1979) was stopped before the self-similar 
structure developed well. 

We dehne an anisotropy parameter A by 


A ^ 2 - 27 / 7 , 


(9) 


where at and ar are the one-dimensional tangential and radial velocity dispersions, re¬ 
spectively. [This anisotropy parameter A was used by, e.g., Giersz and Spurzem (1994).] 
Figure 3 shows the evolving radial prohle of A. As the self-similar core collapse proceeds, 
slight anisotropy proceeds to the inner region self-similarly. In the power-law region, 
A = 0.16, or a‘l/a’1 = 0.92. From hgure 5 of Cohn (1985) we hnd the corresponding value 
of A 0.15. Thus Cohn’s calculation is consistent with our calculation in this respect. 
However, there are differences between the two calculations with respect to the develop¬ 
ment of anisotropy in the outer region of the cluster. For example, at the half-mass radius 
and the 90%-mass radius, A=0.44, 1.62 in our calculation, and A=0.28, 1.24 in Cohn’s 
(1985) calculation, at late times of the core collapse. Thus our calculation indicates the 
development of somewhat stronger anisotropy in the halo. 

It is also useful to compare our results with the results of recent anisotropic gaseous 
models. Louis and Spurzem (1991) obtained pre-collapse (and post-collapse) self-similar 
solutions of their anisotropic gaseous models. Moreover, time-dependent pre-collapse 
solutions of the anisotropic gaseous models were shown in hgures 2(a) and (b) of Ciersz and 
Spurzem (1994); these hgures can be compared with our hgures 2 and 3. The anisotropic 
gaseous models include two numerical constants, Aa and A in Ciersz and Spurzem’s (1994) 
notation. The constants Aa and A are related to the time-scales of collisional decay of 
anisotropy and the heat transport, respectively. These are free parameters in the context 
of the gaseous models, and their values can be determined only through comparison with 
other models such as A^-body models or FP models. The values of a and A in the gaseous 
models depend on AaA. Ciersz and Spurzem (1994) plot the values of a and A in the inner 
power-law region as a function of AaA in their hgures 3(a) and (b). In these hgures, Louis 
and Spurzem’s (1991) self-similar solutions are also plotted, and the self-similar solutions 
and the time-dependent solutions coincide rather well. Ciersz and Spurzem (1994) chose 
Aa = 0.1 and A = 0.4977 as standard values through comparison with A^-body and 
isotropic FP models. With these standard values, the gaseous models give a = 2.22 and 
A = 0.26. This value of a is very close to our value a = 2.23, and this value of A is a 
little higher than but not far from our value A = 0.16. 

Figure 4 presents the evolution of the core collapse rate. 


e = c(0) 


(ilnp(O) 

dint 


as a function of the scaled escape energy. 


a;o = 30(O)/n^(O), 


( 10 ) 


( 11 ) 


where tr(0) is the central relaxation time and nm(0) is the central total velocity dispersion. 
Here the central relaxation time fr(0) is dehned as 


c(0) 


0.065<(0) 
G^mp(O) In A’ 


( 12 ) 


where m is the mass of a star and In A is the gravitational Coulomb logarithm (Spitzer, 
Hart 1971a; Spitzer 1987). The parameter A can be expressed as A = fiN, where N is 
the total number of stars in the cluster and p is a numerical constant (Spitzer and Hart 
adopted p = 0.4). In the present problem, we do not have to specify the values of p and 
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N, because these (in fact, term iV/lnA) can be eliminated from the FP equation if we 
use the time measured in units of the relaxation time. (However, if we include binary 
heating terms into the FP equation, we have to specify the values of /i and N.) The value 
of Xq monotonically increases as the core collapse proceeds. Cohn (1980) found that Xq 
tended to a limiting value of about 13.9, while we hnd that xq ~ 13.5 at the end of our 
calculation. If we continued the calculation further, xq will increase beyond 13.5. Our 
calculation gives the asymptotic constant value of .^ = 2.9 x 10-3, while Cohn’s (1979) 
2D calculation gave ^ = 6.0 x lO-^ and Cohn’s (1980) ID calculation did ^ = 3.6 x lO-^. 
Heggie and Stevenson (1988) obtained the value of .^ = 3.64 x lO-^ for the pre-collapse 
self-similar solution of the ID FP equation. 

While the value of ^ in our anisotropic model is smaller than that in the isotropic 
models, the value of ^ in Cohn’s (1979) anisotropic model is larger. Which tendency 
is true? Spitzer (1987, p.95) argued that the higher value of ^ in anisotropic models 
was also implied by a Monte Carlo calculation by Duncan and Shapiro (1982) Spitzer 
estimated the value of ~ 8 x IQ-^ for that calculation. However, the asymptotic value 
of ^ in Cohn’s (1979) calculation may be less reliable than that in our calculation, because 
Cohn’s calculation was accompanied by large error in energy conservation (~ 11%) and 
because the calculation did not follow the core collapse so deeply. The reliability of the 
value given by the Monte Carlo calculation is further lower. On the other hand, the lower 
value of ^ in anisotropic models was argued by Louis (1990). Louis obtained the self¬ 
similar solutions of isotropic and anisotropic fluid-dynamical models for star clusters and 
found the smaller value of ^ in the anisotropic model. He discussed the physical reason 
for the decrease of the core collapse rate in anisotropic clusters. He derived a relation 
among three different types of deviations from thermal equilibrium, namely the dehciency 
of high-velocity stars, the temperature gradient and anisotropy: the temperature gradient 
can be balanced by the deficiency of high-velocity stars which causes a non-zero collision 
term and thus core collapse, and, for the anisotropic case, by anisotropy. Louis argued 
that generation of anisotropy may be interpreted as an attempt to avoid core collapse. 
Louis’s discussion may be expressed in another way as follows. In anisotropic (i.e. real) 
clusters, the heat is transferred from the core more by stars on more radial orbits than by 
stars on more circular orbits. In isotropic clusters, however, such a bias to radial orbits 
cannot occur, and high-energy circular orbits are enforced to be produced more than in 
anisotropic clusters. High-energy stars on more radial orbits interact with the more inner 
region than stars on more circular orbits with the same energy, or, roughly speaking, 
the high-energy stars on more radial orbits still (in part) stay in the core. This implies 
that the heat is transferred from the core to the outer region more quickly in isotropic 
clusters. In other words, we may say that the reduction of the 2D FP equation to the 
ID FP equation by averaging over angular momentum space causes artihcial diffusion 
in addition to real diffusion. Consequently the core collapse proceeds faster in isotropic 
clusters. 

Figure 5 shows the evolution of the quantity, 

_ (ilnu^(O) 

(ilnp(O) ’ 

which specihes the equation of state of the core. The small fluctuations in the curve may be 
due to numerical error. Our calculation gives the asymptotic constant value of 7 = 0.10, 


(13) 
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while Cohn’s (1979, 1980) 2D and ID calculations gave 7 = 0.12,0.10, respectively. The 
coincidence of the asymptotic value of 7 = 0.10 between our 2D calculation and Cohn’s 
(1980) ID calculation indicates that anisotropy does not affect the equation of state of the 
core (in fact, exact isotropy is always established inside the core). In a simple homological 
model of core collapse (Lynden-Bell, Eggleton 1980), 7 is related to a as 7 = (a — 2)/a. 
This relation is well satished with the values a = 2.23 and 7 = 0.10. The time remaining 
to complete core collapse, r, is represented by a hxed multiple of the current central 
relaxation time [equation (13) of Cohn (1980)]: 

^ = [(1 - ^7)^]"^fr(0). (14) 

The asymptotic constant values of ^ and 7 in our calculation gives r ~ 410tr(0), while 
Cohn (1980) obtained r ~ 330tr(0). 

Figure 6 a shows the evolution of Lagrangian radii containing inner 1, 2, 5, 10, 20, 30, 
40, 50, 75, and 90% of the cluster mass. The result of a ID FP calculation is also shown 
by dotted curves. The time is measured in units of the initial half-mass relaxation time 
trh,i- The half-mass relaxation time frh is dehned by 


^ 1 / 2 ^ 3/2 


(15) 


where rh is the radius containing half of the total mass (Spitzer, Hart 1971a; Spitzer 
1987). We hud again slower core collapse in the 2D calculation from figure 6 a. The 2D 
calculation gives the core collapse time of fcoii = 17.6trh,i and the ID calculation does 
tcoii = 15.6trh,i- To clarify differences between the 2D and ID models except in the time 
scale, we also present hgure 6 b where the time axis of ID calculation is multiplied by a 
constant factor so that the collapse time in the ID calculation should coincide with that 
in the 2D calculation. In this hgure, we do not hud any signihcant differences between 
the two models for 1-75% Lagrangian radii. However, the 90% radius of the 2D model 
expands further than that of the ID model, that is, the 2D model has more extended 
halo. 

The density prohles at epochs when the central density increases by about 4.5 orders 
of magnitude are shown in hgure 7 for the 2D model (the solid curve) and the ID model 
(the dotted curve). We hud again that the density prohles for ID and 2D models are 
very similar inside the half-mass radius, but that the 2D model has the higher density 
in the outer halo. The density prohle in the outer halo is approximated by a power law 
p oc rather well for the 2D model, while no corresponding power law is observed for 
the ID model. (The initial density prohle in the halo is asymptotically given by p oc r~^ 
as indicated by Plummer’s density law.) This power law is not strictly established, but 
the value of the logarithmic density gradient is between —3 and —4 in the halo (hgure 
2). Spitzer and Shapiro (1972) demonstrated by using the (F, J)-space FP equation that 
the steady state density prohle in the halo is given by p oc Spitzer and Shull 

(1975) showed this power law was observed in their Monte-Carlo calculation started from 
Plummer’s model. However, Cohn (1979) did not observe this power law. He stated this 
discrepancy might be due to the use of an energy cutoh in his calculation, corresponding 
to a tidal limiting sphere of radius rt = lOOrg {tq is a length scale parameter of Plummer’s 
model); the density must vanish at r = rt in his calculation, so that din p/d In r must 
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drop considerably near r = r^. In the present calculation we used different boundary 
conditions: the flux Fg = 0 on an energy cutoff boundary E = Fmin. Therefore the 
density does not have to vanish at the outermost radial grid point. 

The well-developed halo in the 2D model is due to the emergence of stars on very 
eccentric orbits from the core into the halo. The halo is dominated by these eccentric 
orbits. This is clearly seen in hgure 8 which shows the radial and tangential velocity 
dispersion prohles for the 2D model at the same epoch as in hgure 7. In the halo the 
radial velocity dispersion exceeds the tangential velocity dispersion considerably. We hnd 
from hgure 3 that the value of the anisotropy parameter A increases monotonically and 
approaches to the maximum value of 2 as the radius increases. A power law prohle 

oc in the outer halo, which corresponds to the constant mean square angular 
momentum, were indicated by several authors (e.g. Spitzer, Hart 1971b; Henon 1971). 
We hnd from hgure 8 that this power law gives a reasonable ht to the result of the 2D 
FP calculation. (The initial velocity dispersion prohle in the halo is asymptotically given 
by p oc r“h) 

Figure 9 shows the evolution of the mass-averaged velocity dispersions between 75% 
and 90% Lagrangian radii. Giersz and Heggie (1994a) showed the same hgure which 
was produced from their statistical data of 1000-body calculations. The evolution of the 
velocity dispersions in our model seems to be qualitatively similar to that in the 1000-body 
model. However, there are some quantitative diherences: e.g., in the 1000-body model 
~ 0.075 at the collapse time, while in our 2D FP model ~ 0.055, 2al ~ 0.04 
at the collapse time. Therefore the value of the anisotropy parameter A is higher in the 
2D model. We note the core collapse was not actually completed in the 1000-body model; 
the core contraction stopped due to binary heating when the central density increased 
only by about two orders of magnitude. The lower value of A in the 1000-body model 
may be due to this incomplete core collapse. 

In hgure 10a, the distribution function / at the same epoch as in hgure 7 is plotted as 
a function of the scaled angular momentum R. Diherent curves correspond to diherent 
energies (see the hgure caption): the upper curves correspond to larger E (lower energy 
in a usual sense) and the lower curves to smaller E (higher energy). For large F, the 
distribution function / is actually independent of the angular momentum F, and is very 
close to an isotropic Maxwellian one. For small F, / increases very rapidly as R goes to 
zero, and / for larger R has been almost unchanged from its initial value. As described 
above, the halo is predominated by stars of very eccentric orbits (F ~ 0). These stars 
come from the inner region of the cluster. This is evident from hgure 10b where / is 
plotted as in hgure 10a, but the abscissa is the pericenter radius rp(F,F). For small F, 
/ has large values when rp lies within the inner region of the cluster (the initial half-mass 
radius is rh,i ~ 0.77). 

In hgures 6 and 9, the time is measured in units of the initial half-mass relaxation 
time. To obtain the physical time, we must specify the number of stars N and the value of 
the coefficient fi in the Coulomb logarithm ln{fiN). While a conservative value is /i = 0.4 
(Spitzer, Hart 1971a), Giersz and Heggie (1994a) obtained /r = 0.11 by comparing A^-body 
models with isotropic gaseous and isotropic FP models. Various estimates of /i by other 
authors are listed in table 2 of Giersz and Heggie (1994a). Rough comparison between our 
2D FP calculation and Giersz and Heggie’s (1994a) 1000-body calculation with respect to 
the evolution of the Lagrangian radii indicates that the value of /r = 0.11 gives a reasonably 
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good agreement of the two results. Since the core collapse proceeds more slowly in the 2D 
model than in the ID model as described above, the value of /i somewhat larger than 0.11 
may give a better agreement. However, comparison with A^-body models is complicated by 
the inevitable formation of binaries in iV-body models; the core collapse stops due to the 
binary heating at relatively early stages in small A^-body systems. Detailed comparison 
as Giersz and Heggie (1994a, b) did is necessary to determine the best empirical value of 
yU. The initial half-mass relaxation time is frh,i = 19.67 for /i = 0.11 and N = 1000. 

6 Summary 


We have developed the numerical codes to solve the orbit-averaged {E, J)-space FP 
equation with high accuracy and have investigated in detail the pre-collapse evolution 
of single-mass spherical star clusters where velocity anisotropy is allowed. The purpose 
of this work is to improve Cohn’s (1979) work in order to understand the evolution of 
anisotropic star clusters thoroughly. Anisotropic (2D) FP models have been seldom used 
for the last hfteen years because of the difficulties of the numerical calculations. In 
particular the accuracy of total energy conservation was bad in 2D FP calculations. It 
was thought that the origin of this error was in the integration scheme of the 2D FP 
equation. 

In this paper we have developed two different integration schemes. One is the hnite 
difference scheme where the Chang-Cooper scheme is simply applied only for the energy 
direction. The other is the hnite element scheme where the test and weight functions 
suggested by the generalized variational principle are used. Using these schemes we could 
follow the core collapse until the central density increased by about 14 orders of magnitude. 
The total energy error was about 1% at the ends of these calculations. We believe that the 
accuracy of the present calculations is enough to trust the quantitative results presented 
in section 

The results obtained by the hnite diherence and the hnite element schemes are gen¬ 
erally in good agreement. Although the hnite diherence scheme seemed to be better to 
follow the core evolution at late stages of the core collapse, we could not dehnitely conclude 
which of the two schemes was better through the whole evolutionary stages. The accuracy 
of the hnite element scheme may be improved by using higher-order basis functions. The 
sources of the numerical errors are not only in the integration of the FP equation but 
also in the calculation of the dihusion coefficients and in the potential-recalculation steps. 
In these calculations many multi-dimensional numerical integrations and interpolations 
are involved, which are of course more difficult to perform than one-dimensional ones. If 
these calculation schemes are improved, the numerical errors will be reduced further. 

Self-similar evolution of the core collapse is clearly seen in the present calculations. In 
the self-similar region between the isothermal core and the outer halo, the density prohle 
is characterized by a power law p oc and velocity anisotropy is characterized by a 

constant ratio oljal = 0.92. This power-law index of the density prohle is the same as 
that in the isotropic FP model. As the core collapse proceeds, the core radius and mass 
go to zero and the velocity anisotropy penetrates into the inner region. Therefore, when 
the core collapse completes, the velocity distribution is anisotropic everywhere in the 
cluster. At late stages of the core collapse, the collapse rate ^ = C(0)(ilnp(0)/(ilnf and 
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the quantity 7 = (ilnu^(0)/(ilnp(0) tend to asymptotic constant values .^ = 2.9 x 10“^ 
and 7 = 0.10. The asymptotic values in the isotropic model are ^ = 3.6 x 10“^ and 
7 = 0.10. With these values of ^ and 7 , we can estimate that the time remaining until 
complete collapse, r, is r = 410tr(0) in the anisotropic model, and r = 330tr(0) in the 
isotropic model. The collapse time is 17.6trh,i in the anisotropic model and 15.6trh,i in the 
isotropic model, when Plummer’s model is chosen as the initial condition. 

In the outer half of the mass, strong velocity anisotropy has developed. At late times 
the ratio of the tangential velocity dispersion to the radial one decreases monotonically to 
nearly zero as the radius increases: e.g., cr^/a^ ~ 0.78, 0.19 at the half-mass and the 90%- 
mass radii, respectively, at the end of the collapse. The tangential velocity dispersion 
prohle in the outer halo is reasonably approximated by a power law oc r“^. This 
indicates that the mean square angular momentum is almost constant independent of the 
radius. Although the density prohles for ID and 2D models are very similar inside the 
half-mass radius, the 2D model has the higher density in the outer halo. The density 
prohle in the outer halo is approximated by a power law p oc which was indicated 

by Spitzer and Shapiro (1972). The well-developed halo in the 2D model is dominated by 
eccentric orbits. These eccentric orbits have been produced by strong relaxation in the 
core. 

In this paper the detailed direct 2D FP calculations have been performed. We have 
seen that the relaxation process is always accompanied by the velocity anisotropy produc¬ 
tion. We have found that the core collapse in single-mass clusters is not seriously affected 
by anisotropy, though the core collapse rate is reduced a little by the development of 
the anisotropy. However, we do not know very much the effects of anisotropy on the 
post-collapse evolution, and on the evolution of more realistic clusters (e.g., multi-mass 
clusters, tidally limited clusters, etc.). The studies of the effects by using anisotropic 
gaseous models have been progressing (Giersz, Spurzem 1994; Spurzem 1994; Spurzem, 
Takahashi 1994). On the other hand, the 2D FP models have been seldom used during 
the last hfteen years, mainly due to the numerical difficulties. The present paper has 
shown that 2D FP calculations can be performed with reasonable numerical accuracy. 
We believe that we do not have to hesitate to do 2D FP calculations now. The evolution 
of more realistic cluster models will be studied in future works. 
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differencing scheme. I also thank Dr. M. Itoh for giving me his 2D FP solver which is a 
original of the present hnite-difference scheme. The early part of this work was carried 
out during my stay at the Institute for Theoretical Physics, Santa Barbara, California, 
USA, supported by the NSF under Grant No. PHY89-04035. This work was supported 
in part by the Grand-in-Aid for Encouragement of Young Scientists by the Ministry of 
Education, Science and Culture of Japan (No. 1338). 

Appendix 1 Discretization Schemes 
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In this appendix we describe the discretization schemes of the FP equation (^) in 
detail. Two independent variables are denoted by (x, y) here. The discretization is done 
by using the hnite difference method and the hnite element method. 

Appendix 1. 1 The Finite Difference Scheme 


We represent f{x,y,t) by the discrete set = f{xi,yj,tn)- We write Ax, Ay, and 
At for X, y, and time mesh sizes, respectively: Xj+i = Xi + Ax, yj+i = yj + Ay, and 
tn+i =tn + At. We do the hnite differencing of the FP equation as 


An 


_ Jtn 

J i,i Ji.j 








F ■ ■ 1 


At 


Ax 


Ay 


(Al) 


where 



F ■ 1 ■ = 

Xl+^,J 


-D^ .^1 . 

XXl+^,J 

F ■ ■ 1 = 

_ pin £ 

^yid+y id+^ 

-D^ . .^1 
y^^,J+2 


fi+l,j fi,j _ jjn 

Ax 

+ | _ jjn ~ fi,j 

Ax Ay 


(A2) 

(A3) 


and similar expressions for the other hux terms. Here, the diffusion coefficients are evalu¬ 
ated at t = tn, and fj is some linear combination of /"■ and In this paper, we adopt 

the Crank-Nicolson scheme: fj = [ff + ff^^)/2. Hereafter we omit the tildes. In the 
right side of equation (|A2|) , j_,_i is given by 

/i+i j+i = ^(/ij + fi+i,j + + /i+i,i+i), (A4) 


and j_i is given by a similar expression. We represent in the hrst term of 


*+o J 


equation (|A2|) as 


A (1 ^xi,j) 


(AS) 


As in the Chang-Cooper scheme, a key point is a choice of the value of parameter 
We may simply choose = 1/2. However, imitating the Chang-Cooper scheme, we 
determine its value as follows. First we consider the equilibrium solution of the FP 
equation. If we assume that / is independent of y and that the diffusion coefficients are 
independent of /, the solution of = 0 is 


f = Cexp 




-dx 


(A6) 


where C is a constant of integration. Then, in the equilibrium state. 


/^+ij 

fid 


= exp 


exp 


r^i+i D^{x',yj,t) 
JXi Dxx , yj, t) 

D 


dx 


xi+yd 


^xxi+^J 


Ax 


(A7) 
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We choose 6xij so that this relation will be reproduced by setting the hnite-differenced 


form of the flux 


to be zero with the neglect of the third term. Then we hud 


^xi,i 


1 


Wxij exp{wxi,j) - I' 


(A8) 


where 


Ax 


^xxi+^,j 


(A9) 


The parameter 6 monotonically decreases from 1 to 0 as tc goes from — cx) to +cxo, and 
6 = 0.5 for tc = 0. Thus the value of 6 is signihcantly different from 0.5 where the 
advection dominates over the diffusion, i.e. \wx\ is large. We can take the same procedure 
for the hnite-differencing of Fy. 

This scheme is not a real 2D generalization of the Chang-Cooper scheme. In fact, it 
does not ensure non-negativeness of the distribution function, nor does it preserve the 
general equilibrium solution. However, if the problem is reduced to the ID problem (i.e. 
/ depends only on x or y) in some region, the present scheme is virtually reduced to 
the (ID) Chang-Cooper scheme. In the present problem, since an isotropic Maxwellian 
distribution is almost established in the core, our scheme is expected to work well there. 
In fact the adoption of our scheme reduced the total energy error by a factor of more than 
ten relative to the simple centered differencing. Numerical experiments revealed that this 
Chang-Cooper-like differencing for i?-direction did not affect the results, and that Sr was 
always close to 0.5 everywhere. Therefore we simply set 5r = 0.5 in the standard runs. 

Cohn (1985) also reported that he investigated several alternative generalizations of 
the Chang-Cooper scheme and that all of these improved energy conservation. Although 
his schemes were not explained in detail, our scheme is possibly similar to them. 


Appendix 1. 2 The Finite Element Scheme 


We approximate the solution / by using a set of basis functions {ipi, ip 2 ^ ■ ■ ■ i^n} as 

N 

g{x, y, t) = log /(x, y,t) = Y, 9j{t)(Pj{x, y) ■ (AlO) 

i=i 

Applying the method of weighted residuals, we multiply equation (^) by weight ipi and 
integrate it over the entire domain D. Thus we obtain 

IL li 

The right side of this equation can be rewritten by application of Green’s theorem, which 
produces 

= + (A12) 

where n denotes the unit normal to the boundary T of D, directed outward. The second 
term of the right side of this equation vanishes if (pj = 0 or F ■ n = 0 on the boundary 
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r. Because this condition holds true in the present case, we omit this term hereafter. 
Substituting equation (|A10|) into equation (|A12|) , we obtain 


= I] + k , (A13) 

i 3 


where 


rriij = 


h-■ — 


h = 




difj 


dx 


dipj\ d(p. 


D^rr ^ ^ + Dxv n ^ 1 -1" I D 


dx 


fdxdy. 

in \ ox oy j 


d^j , _ dipj\ dipi 


yx 


dx 


+ D 


yy 


dy J dy 


f dx dy, 

(A14) 


Then equation (A13|) is discretized with respect to time as 


E 


m 


n _ 

ij 


^3 


At 


= E *^5® +. 


(A15) 


where gj is a linear combination of and g]'~^^- In fact we take gj = {gj'~^^ + g])/‘^- We 
simply evaluate rriij, hj, and k at t = tn- 

We use rectangular elements to discretize the domain ff . Node points are located 
at the four corners of the element. Two-dimensional piecewise bilinear polynomials are 
used as the basis functions. Concerning the more detailed procedure of the hnite element 
method, see appropriate textbooks (e.g., Fletcher 1991). 

We may expand solution / as / = Y. instead of equation (|A10|) . Test 

calculations by using the hnite element scheme with this trial function were also carried 
out. The mass and energy errors were larger than in the calculation with the trial function 
( |A1(J|) . The value of / varies by many orders of magnitude over the domain where the 
equation is solved. This indicates that we had better use In / as a variable. This may 
be one simple reason why the trial function ( |A1U|) is better. We also note that the form 
of the trial function (|A10| ) is naturally derived from Inagaki and Lynden-Bell’s (1990) 
generalized variational principle. Taking the variation of the local potential <F [equation 
(5.5) of Inagaki and Lynden-Bell (1990)], we obtain 


= JJ^A^dlnfdxdy 

= 0 . 



+ Fy^^^] dxdy 


dx 


dy 


(A16) 


When we adopt the trial function (|A10|) , we obtain equation (|A12|) again. 
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Figure Captions 

Fig. 1. Evolution of the density profile. The evolving density prohle is plotted every 
200 potential-recalculation time steps. The central density increases and the core radius 
decreases with time. 

Fig. 2. Evolution of the radial profile of the logarithmic density gradient d In p/d In r = 
—a. The prohles at the same epochs as in hgure 1 are plotted. The power-law region 
with a = 2.23 extends into the inner region as the core collapse proceeds. 

Fig. 3. Evolution of the radial prohle of the anisotropy parameter A = 2 — 2(t//(T^. 
The prohles at the same epochs as in hgure 1 are plotted. The self-similar region with 
A = 0.16 extends into the inner region as the core collapse proceeds. 

Fig. 4. Evolution of the core collapse rate ^ = fr(0)[dlnp(0)/dlnt] as a function 
of the scaled escape energy xq = 30(0)/n^(0). The asymptotic constant value of ^ is 
^ = 2.9 X 10-T 

Fig. 5. Evolution of the quantity 7 = dlnn^(0)/dlnp(0) as a function of xq. The 
asymptotic constant value of 7 is 7 = 0.10. The small huctuations in the curve may be 
due to numerical error. 

Fig. 6 . (a) Evolution of Lagrangian radii containing inner 1, 2, 5, 10, 20, 30, 40, 50, 
75, and 90% of the cluster mass. The solid curves are the result of the 2D FP calculation, 
while the dotted curves are that of ID calculation. The time is measured in units of the 
initial half-mass relaxation time trh,i- (b) Same as (a), but the time axis of ID calculation 
is multiplied by a constant factor so that the collapse time in the ID calculation should 
coincide with that in the 2D calculation. 

Fig. 7. The density prohle at the 400th potential-recalculation time step (t = 17.5trh,i) 
in the 2D FP calculation is shown by the solid curve. The dotted curve is the density 
prohle in the ID calculation at an epoch when the central density is about the same as in 
the 2D calculation (t = 15.6trh,i)- The asymptotic line p oc r~^'^ is shown for comparison. 

Fig. 8 . The velocity dispersion prohles at the same epoch as in hgure 7. The solid 
and dashed curves are the radial and ID tangential velocity dispersions in the 2D FP 
calculation, respectively. The dotted curve is the ID velocity dispersion prohle in the ID 
FP calculation. The asymptotic line oc is shown for comparison. 

Fig. 9. Evolution of the velocity dispersions which are mass-averaged between 75% 
and 90% Lagrangian radii. The solid and dashed curves are the radial and 2D tangential 
velocity dispersions, respectively. 

Fig. 10. The distribution function / at the same epoch as in hgure 7 is plotted 
as a function of (a) the scaled angular momentum R, and (b) the pericenter radius rp, 
for a hxed energy E. Diherent curves correspond to diherent energies 77=0.0170, 0.0380, 
0.0847, 0.188, 0.410, 0.869, 1.74, 3.14, 4.90, and 6.54, from the bottom to the top. The 
central potential is 0(0) = 7.69. 
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